This invention relates to estimating the probability of future values of a substantially periodic phenomenon. More specifically, the invention relates to an apparatus and method for making such estimations in ordering the views of a nuclear magnetic resonance (NMR) scan to reduce motion artifacts caused by periodic physiological motion.
NMR has been developed as an imaging method useful in diagnostic medicine. In NMR imaging, as is understood by those skilled in the art, a body being imaged is held within a uniform magnetic field oriented along a z axis of a Cartesian coordinate system.
The net magnetizations of the nuclei of the body are excited into precession by means of a radio frequency (RF) pulse and the decaying precession of the spins produces an NMR signal. The amplitude of the NMR signal is dependant, among other factors, on the number of precessing nuclei per volume within the imaged body termed the "spin density".
Magnetic gradient fields G.sub.x, G.sub.y, and G.sub.z are applied along the x, y and z axis to impress position information onto the NMR signals through phase and frequency encoding. A set of NMR signals may then be "reconstructed" to produce an image. Each set of NMR signals is comprised of many "views", a view being defined as one or more NMR signal acquisitions made under the same x and y gradient fields.
The length of time needed to acquire a sufficient number of views to reconstruct an image may be on the order of several minutes. Therefore, motion of the subject, including cardiac and respiratory motion, is inevitable. Such motion may produce blurring of the image or other image artifacts depending on the image reconstruction technique used.
One NMR image reconstruction technique, associated with "spin warp" imaging, will be described herein. It should be recognized however that the present invention may be advantageously practiced with other NMR image reconstruction techniques, such as multiple angle projection reconstruction as disclosed in U.S. Pat. No. 4,471,306.
Referring to FIG. 1, a typical "spin echo" pulse sequence for acquiring data under the spin warp technique includes: 1) a z-axis gradient G.sub.z activated during a first 90.degree. RF pulse to select the image slice in the z axis, 2) a y-axis gradient field G.sub.y to phase encode the precessing nuclear spins in the y direction, and 3) an x-axis gradient G.sub.x activated during the acquisition of the NMR signal to frequency encode the precessing nuclear spins in the x direction. Two such NMR acquisitions, S.sub.1 and S.sub.1 ', the latter inverted and summed with the first, comprise the NMR signal of a single view "A" under this sequence. Note that only the y gradient field G.sub.y changes between view "A" and subsequent view "B". This pulse sequence is described in detail in U.S. Pat. No. 4,443,760, entitled: "Use of Phase Alternated RF Pulses to Eliminate Effects of Spurious Free Induction Decay Caused by Imperfect 180 Degree RF Pulses in NMR Imaging", and issued Apr. 17, 1984 and assigned to the same assignee as the present invention
Referring to FIG. 2, the acquisition of views to construct an image under spin warp imaging is illustrated by referring to the concept of k-space whose coordinates are k.sub.x and k.sub.y. A two-dimensional slice of the imaged object 10 is subject to the pulse sequence of FIG. 1. The NMR signal comprising one view provides the values within a bounded area of k-space 12 along one "row" 14 of k-space 12 orientated along the k.sub.x axis. The particular row 14 is determined by the relative amplitude of G.sub.y for that NMR signal acquisition: for example, more positive y gradients (each measured by area) will produce the NMR signals filling higher row "numbers".
Under the spin warp imaging technique, once all the rows of the entire k-space area 12 are filled by repeated NMR signal acquisitions at different gradients G.sub.y, the k-space area 12 is subject to a two-dimensional inverse Fourier transform along the k.sub.x and k.sub.y axes to produce the image 16.
The data acquired during a single view may be considered essentially instantaneous compared to the periods of cardiac and respiratory motion for typical NMR studies. The data for a given column in k-space, however, is acquired over many views and the time between the acquisition of the data at the top of a column in k-space and the bottom of the column may be substantial. For this reason, motion artifacts in spin warp imaging are generally disposed along the y-axis of the image corresponding to motion induced variations in the columnar k-space data
Referring still to FIG. 2, if the imaged object 10 is stationary, the data read along a column of k-space 18 has frequency components dependant principally on the distribution of spin densities of the imaged object 10 in the y direction and the change in G.sub.y between views. If the imaged object 10 is moving, however, the data in each column 18 is modulated by the frequency of the movement and hence has frequency components corresponding to the modulation sidebands. These sideband frequencies will be discussed below.
The modulation of the data in a column of k-space by movement of the imaged object can best be understood by considering a single area element "a" in the slice of imaged object 10. As object 10 moves, new material will move into area "a" possibly with a higher or lower spin density. This change in spin density in turn will amplitude modulate the data in k-space 12. Of course, if the areas around area "a" in the imaged object 10 have approximately the same spin density as area "a", no modulation will occur. This corresponds with the observation that motion artifacts are most concentrated around areas of abrupt changes in spin density. Additionally, as object 10 moves, signal from area "a" may have a different spin phase due to, for example, motion in the presence of field gradients. This change in spin phase will phase modulate the data in k-space.
Both amplitude modulation and phase modulation of the data in a column of k-space 18 by motion of the imaged object 10 create frequency side bands understood by those familiar with the principles of modulation. These side bands, when transformed by the image reconstructing Fourier transforms to produce image 16, become visible as repetitions or "ghosts" on either side of high contrast image elements "a". Generally, as discussed, these ghosts are displaced in the y direction.
U.S. Pat. Nos. 4,663,591 and 4,706,026 both entitled: "Method for Reducing Image Artifacts due to Periodic Signal Variations in NMR Imaging", issued on May 5, 1987 and Nov. 10, 1987, respectively, and incorporated herein by reference, describe ghost artifacts resulting from periodic motion in the image object and disclose a method of reducing such artifacts by changing the order in which NMR views are acquired.
The order in which NMR views are acquired and hence the order in which the rows of k-space are filled depends on the order in which the gradient amplitudes of G.sub.y are applied during the NMR signal acquisition. By controlling the order of filling the rows in k-space 12, with regard to the modulation of the NMR signal by the motion of the imaged object 10, the apparent modulation of a column of k-space data 18 may be radically altered.
Specifically, if the k-space rows 14 are acquired so that the row number corresponds with the degree of modulation of the NMR signal, the apparent modulation of the data in a k-space column 18 will be decreased in frequency. This ordering strategy is termed a "low sort" and results in a compression of the ghost artifacts toward area "a". This compression will allow distant structures to be more easily viewed.
Conversely, if the k-space rows are acquired so that odd rows correspond to large values of modulation and even rows correspond to low values of modulation the apparent modulation of the data in a k-space column 18 will be increased in frequency. This ordering strategy is termed a "high sort" and results in an expansion of the ghost artifacts away from area "a". If this expansion is sufficient, it may move the ghosts out of the field of view of the image 16, effectively causing the ghosts to disappear.
It should be noted that the sorting order of the rows: "high sort" or "low sort", must be determined at the time of the NMR data acquisition. The rows may not be collected and then rearranged subsequently for the best sort because G.sub.y, which defines the row number, is necessarily fixed at the time the NMR data is acquired. It is necessary, therefore to anticipate the phase of the modulating physiological motion in advance to determine which row to acquire.
U.S. Pat. No. 4,720,678 entitled: "Apparatus and Method for Evenly Distributing Events over A Periodic Phenomenon" issued Jan. 19, 1988, and incorporated herein by reference, describes several methods of predicting the phase of a generally periodic physiological motion. These methods are adapted for use in selecting row ordering in an NMR scan so as to produce "high" or "low" sorting as described above. One such method, as applied to respiration will be briefly summarized as follows:
Referring to FIG. 3(a) a representative respiration waveform 20, as may be produced by a belt type respiration transducer measuring chest wall extension, has a periodic amplitude variation y(t). This waveform 20 is sampled for a period of time to determine its extreme values and the probability density function of its values over time. The probability density function of waveform 20 is simply a measure of the relative probability that the respiration waveform 20 will have a value within a given amplitude range. For the generally sinusoidal respiration waveform 20 shown in FIG. 3(a), the probability density function will be highest for values near the waveform's crest and trough because a relatively greater amount of time is spent by the waveform 20 at these values than at values at the center of the waveform 20. In practice, a probability density function may be acquired by tallying the particular values of a waveform on a histogram whose horizontal axis defines intervals of waveform value and whose vertical axis defines number of occurrences divided by total occurrences.
The probability density histogram of respiration waveform 20 may be integrated along its horizontal axis to produce an "integral histogram" which indicates the probability that the waveform 20 will have a value equal to or less than the value indicated on the horizontal axis. For purposes of graphical clarity, the integral histogram 22, shown in FIG. 3(b), has been rotated so that its vertical axis is waveform value, to coincide with the vertical axis of FIG. 3(a) which also shows waveform value. The horizontal axis of the integral histogram 22 of FIG. 3(b) shows probability.
The integral histogram 22 serves as a transfer function that maps values of the respiration waveform 20, whose values vary in probability of occurrence, to a triangle waveform 24, shown in FIG. 3(d), whose values are equally probable. The triangle waveform 24 will be termed phase and denoted .phi.(y) because it is a linear function of time related to the amplitude of the respiratory waveform y(t). The phase .phi.(y) is used to determine which view or row of k-space should be next acquired so as to reduce any motion artifacts as will be described further below.
Referring to FIG. 3(b), prior to the NMR signal acquisitions, each value G.sub.y (i) is assigned to an evenly spaced gradient reference phase value .phi..sub.ref (i) shown along the horizontal axis of the integral histogram 22. For the low sort to be illustrated, negative gradient values G.sub.y (i) are assigned to low phase values .phi..sub.ref (i). Hence, in this example, i may be considered a row number in k space. For clarity only eleven reference phase and gradient values are shown.
Immediately before an NMR signal acquisition, the value of the respiration waveform 20 is sampled, the closest reference phase value .phi..sub.ref (i) is determined from the histogram 22, and the associated G.sub.y (i) value is chosen. Referring to FIG. 3(a) at time t.sub.1 the respiration waveform 20 is at an amplitude y(t.sub.1) indicated by point "A". Referring to the integral histogram 22 of FIG. 3(b) the amplitude of sample "A" corresponds to reference phase .phi..sub.ref (0) and the gradient value G.sub.y (-5). The NMR signal is thus acquired with gradient G.sub.y (-5) and stored in the proper elements of an array representing the zero-th k-space row 14.
Immediately before the time of the next NMR acquisition, t.sub.2, the respiration waveform 20 is sampled at point "B" whose amplitude corresponds to a reference phase .phi..sub.ref (3) and a gradient G.sub.y (-2). Again the NMR signal is acquired and the third k-space row 14 filled. This process is continued until each gradient G.sub.y (i) has been used which indicates that the array in k-space 12 is complete. The k-space rows are thus not acquired in numerical order but are "reordered" by the use of the integral histogram 22. The reordered k-space data 12 is ultimately processed by a two dimensional Fourier transform to reveal an image 16 as has been discussed.
The amplitude of the respiration waveform 20 at the time of the acquisition of the NMR signal, and thus for each row of k-space, is shown in the plot in FIG. 3(c). The change in amplitude of the respiration waveform 20 with k-space row number reflects the possible modulation of the data in a k-space column, and hence the plot of FIG. 3(c) is named "modulation in k-space".
As may be seen, acquiring the k-space rows according to the integral histogram lowers the modulation frequency of plot 22 in comparison to that which would be obtained if the k-space rows were acquired in numerical order. As discussed above, this lowering of the frequency of respiration induced modulation will tend to compress the ghost artifacts in the resulting image 16.
This example of the ghost reduction technique of U.S. Pat. No. 4,720,678, cited above, incorporates a number of simplifications. First, the amplitude y(t.sub.i) of the respiration waveform 20, at a given sampling time t.sub.i, may not produce a phase value .phi.(y.sub.i), when transformed by integral histogram 22, that is equal to any reference phase .phi..sub.ref (i). The acquisition of the NMR signal at times t.sub.i occurs at an essentially constant rate, dictated by the particular pulse sequence, and is independent of the amplitude of the respiration waveform 20. In cases where the phase of the sampled value .phi.(y.sub.i) does not match any reference phase .phi..sub.ref (i), the closest reference phase .phi..sub.ref (i) will be chosen.
The second simplification is that the above example ignores the possibility that two samples will have values y(t.sub.i) corresponding to the same reference phase .phi..sub.ref (i). The even spacing of the reference phases values .phi..sub.ref (i) within the range of the equiprobable phase function .phi.(y) reduces the chance that the phase .phi..sub.ref (i) will be selected twice prior to all other phase values .phi..sub.ref being selected once. This is the advantage of the integral histogram transformation. However, if a sampled value of the respiration waveform 20 matches a given reference phase .phi..sub.ref (i) that has already been used, the next closest unused reference phase .phi..sub.ref (i) will be chosen instead.
It may be noted that these problems of matching the values of the respiration waveform y(t.sub.i) to a particular reference phase value .phi..sub.ref (i) arise only because the NMR image must be acquired in a short period of time. If time were no object, the precise value of y(t.sub.i) corresponding to a particular reference phase .phi..sub.ref (i) could be obtained. Thus errors resulting from imperfect matching of the values of the respiration waveform y(t.sub.i) to a particular reference phase value .phi..sub.ref (i) will be termed "sampling limitation errors".
A third simplification in the example of FIGS. 3(a)-(d) is in the uniformity of the respiration waveform 20. Typically such waveforms vary in amplitude, baseline, shape, and frequency. The predictive accuracy of the integral histogram 22 is not affected by changes in respiration waveform frequency but suffers if the amplitude of the respiration waveform 20 as measured from crest to trough, or if the baseline of respiration waveform 20, measured along the troughs or the shape of the waveform, varies during the acquisition of the NMR signals.
Referring to FIG. 4(a), a varying respiration waveform 28 has an amplitude which decreases after the compilation of the integral histogram 22 shown in FIG. 4(b). Applying a "low sort" to the sampled data points in respiration waveform 28, as described above, yields the jagged plot of modulation in k-space 30 of FIG. 4(c). The fluctuations in respiration amplitude values in this modulation plot 30 impairs the effect of the low sort by introducing high frequency modulation into the k-space data.
The jagged peaks in plot 30 result from failure of the histogram 22 to accurately predict the probability of the values of respiration waveform 28 after its drop in amplitude. The failure of prediction results in early gradients G.sub.y (i) being mis-assigned, which in turn, forces later gradients G.sub.y (i) to be assigned to inappropriate respiration values. It will be apparent from this discussion that a similar problem will arise for changes in the baseline of the respiration waveform 28 and that similar problems will also affect "high sort" strategies. The errors resulting from failure of the histogram 22 to reflect the probability of values of the respiration waveform 28 as it changes will be termed "prediction errors" to distinguish them from "sampling limitation errors" described above and from a third type of error to be discussed below.
The prior art has addressed "prediction errors" resulting from unstable respiration waveforms in two ways. First the waveform may be preprocessed by automatic gain control ("AGC") and automatic offset control ("AOC") circuitry so that the integral histogram transfer function is presented with an essentially constant amplitude respiration waveform with constant baseline. Second, the histogram may be continually updated by adding new data to the histogram and discarding the oldest data in a "rolling" fashion. Changes in amplitude or baseline will thus be rapidly accommodated by the integral histogram which changes shape.
Inevitably, however, changes in the integral histogram will change the functional relationship between y(t) and .phi.. Referring to FIG. 5(a), a respiration signal 29, similar to that of respiration waveform 28 of FIG. 4(a), changes in amplitude at time t.sub.c. Two histogram curves 30 and 32 represent the probabilities of the values of the respiration waveform before time t.sub.c and after time t.sub.c respectively. By using histogram 30 for sampled points of respiration waveform prior to t.sub.c and using histogram 32 for sampled points subsequent to time t.sub.c, the goal of a rapidly adapting histogram may be achieved.
FIG. 5(c) shows the modulation of data in k-space based on the view ordering using histograms 30 and 32. Clearly, the matching of respiration waveform values to rows in k-space has significant errors despite the rapid adaptation of the histogram of FIG. 5(b) to the respiration waveform. Although the histogram curves 30 and 32 accurately reflect the probability of values of the respiration waveform on an instantaneous basis, the functional relationship between y(t) and .phi. during the matching of respiration values to gradients G.sub.y is weakened. Points H, I, J, and K map perfectly to reference gradients 1, 4, 7, and 10, on histogram curve 32 and are correctly ordered with respect to each other, however the actual values y(t.sub.i) have no correspondence to the values of y(t.sub.i) selected by histogram 30. These errors in ordering, resulting from the change in the histogram during the matching of values to elements, will be termed "correspondence errors". It is important to note that "correspondence errors" are exacerbated by rapid changes in the histogram curve, the very method that is used to reduce "prediction errors" as described above.
It should be noted additionally that correspondence errors will also occur if the histogram is held fixed and the respiration waveform is preprocessed by AGC and AOC circuitry. AGC and AOC processing accomplishes the same end as does rapid updating of the integral histogram. Both procedures weaken the relationship between y(t) and .phi..